Code
mod_interventions <- c("No intervention", "Naloxone", "Naloxone", "Prescription guidelines", "Safer Supply", "All interventions")
changes <- c("No changes", "")Cost benefit analysis of the federal interventions that aim to reduce the harms of opioid crisis
Final Project | EPIB676
1 Background
—What is the decision problem your model aims to inform? What is already known, and what are the key questions you hope to address? Why is this important?—–
This project aims to evaluate the future economic impact of these interventions by conducting a cost-benefit analysis to examine the difference in healthcare costs between the interventions and their projected effect on reducing opioid-related deaths and other important health outcomes
2 Methods
———–What is your analytic design, setting, outcome(s) estimated, alternatives, population, time horizon? How is your model structured, and what dynamics or quantities are captured? What simplifications have you made to the model? What data are you using to parameterize or calibrate your model (ideally, presented in one or more tables)? What uncertainty analysis are you conducting————–
An open cohort Markov model was used (fig. 1) that consists of 31 states representing the pathways of opioid use in Canada over 15 years (2015 to 2029) with monthly cycles for 15+ Canadian population. It is an open cohort model that takes into account the changes in population over time by adding people to the first state (Pain free, no use) every cycle. The yearly increase was assumed to be constant throughout the year (so each increase was divided by 12, as for years in the future, projections of population size in Canada were used and proportion of those 15+ was estimated based on previous data, then a similar approach was used to estimate the monthly increase in the cohort).
— Talk about the model and probabilities —-
2.1 Alternatives
This model was run for 5 scenarios in regards to the previously mentioned interventions: 1) No interventions, 2) Naloxone only, 3) Prescription guidelines only, 4) Safer Supply only, and 5) all 3 interventions. The effects of each of the interventions were estimated by other team members who conducted a systematic review of the literature (Table … shows a summary of the changes associated with each intervention)
For the naloxone intervention: the probability of transitioning from illicit opioid overdose to OD death was reduced by 5%, and that of transitioning from rx opioid overdose to OD deaths was reduced by 3%. As for prescription guidelines the following changes were applied, transitioning from pain free to acute
2.2 Outcomes
Opioid-related costs from the healthcare perspectives for all states will be used in this analysis. These would include costs for hospitalization, ED visits, physician billing, paramedic services, drugs, and other healthcare system related costs. These were identified from publicly available data, as well as published systematic reviews, studies, and reports. I will also examine epidemiological measures: prevalence of OD, number of OD-related deaths and OD-related brain injury
The primary outcomes are difference in 1) costs, 2) opioid-related overdose deaths (OD deaths) and 3) overall deaths excluding OD deaths between interventions models and no intervention model. The secondary outcomes are weighted yearly prevalence of: 1) severe brain injuries due to opioid overdose, 2) moderate brain injuries due to opioid overdose, 3) substance use treatment and 4) prescription opioid use, in addition to monthly prevalence trends of the 4 outcomes. Primary outcomes were included the probabilistic sensitivity analysis (PSA), while only point estimates were calculated for secondary outcomes.
2.3
3 Results
What are your main findings? Should estimate some outcome(s) for >1 alternative, and potentially conduct incremental analysis (compute ICERS) depending on your framework. Report uncertainty analysis findings to give some idea of the robustness of your findings
4 Expansion Plan
If you were to develop this initial analysis into one fit for submission for peer reviewed publication, what steps would you take? Which aspects of the analysis would need to be improved?
5 Discussion
What are the key takeaways from your initial analysis (and/or, what might the takeaways be if a more complete analysis was conducted)? What is the relationship between your findings and the wider literature or policy landscape? What are the strengths and limitations of your analysis?
Code
library(here)6 Appendix
Code
library(here)
library(RColorBrewer)
source(here("02_scripts/01_fun_data.R"))
calib_target_tbl <- read_excel(here("01_data/markov_model_parameters_preprior.xlsx"),
sheet = "calibration_targets")
theme_set(theme_minimal())6.1 Prevalence of Rx opioids use
Code
# prevalence of people on prescribed opioids
v_yrs_prev_rx <- c(2015:2018)
yrs_prev_rx <- length(v_yrs_prev_rx)
prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("prop_opioids_rx_",
str_sub(v_yrs_prev_rx, 3, 4)))))
tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = yrs_prev_rx,
list(NULL,
paste0("tot_pop_",
str_sub(v_yrs_prev_rx, 3, 4)))))
for (i in 1:yrs_prev_rx) {
yr <- v_yrs_prev_rx[1] + i
prop_opioids_rx_uncalib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_uncalib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(v_yrs_prev_rx,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target, group) %>%
rename(grp = group,
target_val = target) %>%
mutate(year_mon = paste(year, month.abb[7], sep = "_"))
prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Model") %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, target_val) %>%
mutate(grp = "Target"))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
levels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(v_yrs_prev_rx,
each = 12),
month.abb[1:12],
sep = "_"))
p1 <- ggplot() +
geom_point(data = prop_opioids_rx_uncalib_tbl,
aes(x = year_mon, y = target_val,
color = factor(year), fill = factor(year))) +
geom_point(data = prop_opioids_rx_target_tbl,
aes(x = year_mon, y = target_val),
color = "red") +
geom_point(data = prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon),
aes(x = year_mon, y = target_val),
color = "blue")+
xlab("Year_Month") + ylab("Prevalence of prescription opioid use") +
scale_fill_discrete(name = "year")+
scale_color_discrete(name = "year")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
plotly::ggplotly(p1)The red points represent the observed target proportions of Rx opioids from CIHI report, the blue points represent an annual average (weighted by population during that month), the other coloured points represent monthly prevalence predicted from the model
Code
p2 <- ggplot() +
geom_point(data = prop_opioids_rx_uncalib_wei_mean,
aes(x = year, y = target_val, color = factor(grp),
fill = factor(grp))) +
xlab("Year") + ylab("Prevalence of prescription opioid use")+
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# ylim(0,0.2) +
# ggtitle("Comparison between uncalibrated annual weighted average and
# observed targets of prevalence of Rx opioids use")
plotly::ggplotly(p2)6.2 Deaths and OD-Deaths
Code
################ Deaths
# number of deaths
v_yr1_deaths <- c(2016:2020)
yrs_deaths <- length(v_yr1_deaths)
num_deaths_uncalib <- rep(NA, yrs_deaths)
for (i in 1:yrs_deaths){
yr <- (v_yr1_deaths[1] - 1) + i
num_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
}
# number of OD-deaths
v_yr1_oddeaths <- c(2016:2021)
yrs_oddeaths <- length(v_yr1_oddeaths)
num_od_deaths_uncalib <- rep(NA, yrs_oddeaths)
for (i in 1:yrs_oddeaths){
yr <- (v_yr1_oddeaths[1] - 1) + i
num_od_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
deaths_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_deaths", "total_od_deaths")) %>%
mutate(target = ifelse(target == "total_deaths", "Total deaths",
ifelse(target == "total_od_deaths",
"Total opioid-related overdose deaths", NA)),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_deaths_uncalib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
mutate(target = "Total opioid-related overdose deaths")) %>%
mutate(year = c(v_yr1_deaths, v_yr1_oddeaths),
group = "Model"))
p3 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "") +
facet_wrap(~target, scales = "free") +
theme(legend.position="bottom")
plotly::ggplotly(p3)Code
plotly::ggplotly(ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>%
filter(target == "Total deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = ""))Code
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>%
filter(target == "Total opioid-related overdose deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total opioid-related overdose deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = ""))6.3 OAT
Code
# number of oat
v_yr1_oat <- c(2018:2021)
yrs_oat <- length(v_yr1_oat)
num_oat_uncalib <- rep(NA, yrs_oat)
for (i in 1:yrs_oat){
yr <- (v_yr1_oat - 1) + i
num_oat_uncalib[i] <- sum(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}
oat_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% "total_oat") %>%
mutate(target = ifelse(target == "total_oat",
"Total OAT", NA),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
mutate(target = "Total OAT",
year = v_yr1_oat,
group = "Model"))
p3 <- ggplot() +
geom_point(data = oat_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "") +
theme(legend.position="bottom")
plotly::ggplotly(p3)6.4 Trace plot
Code
# build-in color palette
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
c("#084C61", "#DB504A",
"#E3B505", "#4F6D7A",
"#56A3A6"))
names(colors_pal) <- NULL
trace_tbl_inc <- as_tibble(mod_basecase$m_M) %>%
mutate(cycle_num = 0:n_cycles) %>%
pivot_longer(1:n_states, names_to = "state", values_to = "count")
trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_od_deaths
trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_mod_bi
trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_sev_bi
trace_tbl_inc$state <- factor(trace_tbl_inc$state,
levels = v_state_names,
labels = v_state_names)
p4 <- ggplot(trace_tbl_inc, aes(x = cycle_num, y = count, colour = state)) +
geom_line() +
ylab("Count") + xlab("Cycle Number") +
ggrepel::geom_text_repel(data = trace_tbl_inc %>% filter(cycle_num == 180),
aes(x = cycle_num, y = count, label = state),
check_overlap = T, size = 3) +
scale_color_manual(values = colors_pal)
plotly::ggplotly(p4)Code
p5 <- ggplot(trace_tbl_inc %>%
filter(state == "BO_OD_DEATH") %>%
filter(cycle_num != 0) %>%
filter((cycle_num + 1) %% 12 == 1) %>%
mutate(year = c(2015:2029)),
aes(x = year, y = count)) +
geom_line() +
ylab("Total opioid-related overdose deaths") + xlab("Year")
# +
# ggtitle("Total opioid-related overdose deaths over time")
plotly::ggplotly(p5)